Recurrence plots and chaotic motion around Kerr black hole 
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Abstract. We study the motion of charged test particles around a Kerr black hole immersed in the asymptotically uniform 
magnetic field, concluding that off-equatorial stable orbits are allowed in this system. Being interested in dynamical properties 
of these astrophysically relevant orbits we employ rather novel approach based on the analysis of recurrences of the system to 
the vicinity of its previous states. We use recurrence plots (RPs) as a tool to visualize recurrences of the trajectory in the phase 
space. Construction of RPs is simple and straightforward regardless of the dimension of the phase space, which is a major 
advantage of this approach when compared to the "traditional" methods of the numerical analysis of dynamical systems (for 
instance the visual survey of Poincare surfaces of section, evaluation of the Lyapunov spectra etc.). We show that RPs and 
their quantitative measures (obtained from recurrence quantification analysis - RQA) are powerful tools to detect dynamical 
regime of motion (regular vs. chaotic) and precisely locate the transitions between these regimes. 

Keywords: black hole physics, magnetic fields, relativity 
PACS: 97.60.Lf, 98.62Js, 04.70.Bw, 05.45.Pq 



INTRODUCTION 

In this work we continue our effort [1, 2, 3] to understand the dynamic properties of charged test particles being 
exposed to the electromagnetic test fields surrounding compact objects - neutron stars and black holes. As we 
bear astrophysical motivation in our mind we choose such fields and such background geometries which combine 
themselves in reasonable models of real situations occurring in the vicinity of these objects. Survey of the test particle 
trajectories might be regarded as the one particle approximation to the complex dynamics of the astrophysical plasma. 

In this contribution we investigate the motion of charged test particles around the Kerr black hole immersed into the 
asymptotically uniform magnetic field (Wald's solution [4]). This field allows motion in the off-equatorial lobes if the 
parameters of the system are chosen carefully [2] . By off-equatorial lobes we mean closed equipotential surfaces which 
enclose the stable off-equatorial circular orbit (so called halo orbit). They symmetrically intersect poloidal (r, 0) plane 
as seen in Fig. 2. We investigate motion of the charged particles in these lobes. We are particularly curious what is the 
dynamic regime of motion (chaotic versus regular) and how does it change if we alter some of the parameters. Besides 
the standard technique of Poincare surfaces of section we employ recurrence analysis [5] and show that recurrence 
plots might be regarded as an alternative tool to surfaces of section. Quantitative analysis of recurrences appears to be 
very useful tool to locate the transitions between these regimes in the terms of selected control parameter. 



KERR BLACK HOLE IN ASYMPTOTICALLY UNIFORM MAGNETIC FIELD 

The Kerr metric in Boyer-Lindquist coordinates f, r, 6, <j) is given as follows [6]: 

ds 2 = -~ [dt - a sinO d(p} 2 + [(r 2 + a 2 )d(j) - adt} 2 + ^ dr 2 +T.d9 2 , (1) 

where we set 

A=r 2 -2Mr + a 2 , £ = r 2 + a 2 cos 2 0, (2) 

and M stands for the mass of the black hole and a for its specific angular momentum. 

In order to incorporate the large scale magnetic field to our considerations we employ the Wald's test field solution 
[4] . Vector potential of this axisymmetric test field may be expressed in the terms of Kerr metric coefficients as follows: 

At = ^B (g tt p+2ag t ,) -®g tt -@, (3) 



A^ = -B {g^+2ag t ^)-—g,^, (4) 

where Q = Q/M stands for the specific test charge of the black hole (unlike the Kerr-Newman solution it does 
not enter the metric). The terms containing Q in above equations may thus be identified with the components of the 
vector potential of the Kerr-Newman solution [6]. Asymptotic behavior of the components justifies the identification 
of the parameter Bq with the strength of the originally uniform magnetic field into which the Kerr black hole has 
been immersed. Wald [4] has shown that in the case of parallel orientation of the spin and the magnetic field Bo the 
black hole selectively accretes positive charges (negative for the antiparallel orientation) until it is charged to the value 
<2w = 2Boa. Although the use of the Wald's charge gw is not obligatory (since the presence of another charging 
mechanism could be considered for the astrophysical black holes) it may be regarded as a preferred value. 



RECURRENCE PLOTS AND RECURRENCE QUANTIFICATION ANALYSIS 

Besides various "traditional" methods of the numerical analysis of dynamical systems (for instance the visual survey 
of Poincare surfaces of section, evaluation of the Lyapunov spectra [7] etc.) it appears useful to employ rather novel 
approach based on the analysis of recurrences of the system to the vicinity of its previous states. Recurrence plots 
(RPs) as a tool to visualize recurrences of the trajectory in the phase space were introduced by Eckmann et al in 1987 
[8]. RP method is based on the examination of the binary values that are constructed from the phase space trajectory 
x{t). Construction of RPs is simple and straightforward regardless of the dimension of the phase space which is a 
major advantage of this approach. Binary values of the recurrence matrix R, 7 may be formally expressed as follows: 

R ij (e) = ®(e-\\x(i)-x(j)\\), i,j = l,...,N (5) 

where e is a predefined threshold parameter, © represents Heaviside step function and N specifies the sampling 
frequency which is applied to the examined time period of the trajectory x(t ). 

Selection of the norm | . | which should be used to detect recurrences in the phase space is not straightforward. 
Although simple norms like L 2 (Euclidean norm) are usually applied directly in this context, we want to reflect 
the curvature of the spacetime as much as possible. Thus we measure distances in the ZAMO's hypersurfaces of 
simultaneity following standard 3 + 1 splitting procedure [9]. Projection tensor we apply to the phase space constituents 
x 11 , Tip is given as follows: 

7|iV =^ 1 uv+M 1 u«v, (6) 

where g^ v is a spacetime metric (1) and stands for the coordinate components of ZAMO's four-velocity u which 
may be expressed as follows [10]: 

setting A = (r 2 + a 2 ) 2 - a 2 Asin 2 G and Q. = ^Mr. 

Consequent order estimates lead to the conclusion that the computation of the recurrence matrix R, 7 may be carried 

out using Euclidean norm applied to the following set of coordinates (^/gTrf, ^fgeed , ^/g^(p , \fg FT n r , \Z8 de ^e , Tyf^^)' 
provided that the threshold e as the upper limit of the distance at which we shall apply this norm is significantly 
smaller than the characteristic scale of the spacetime curvature. As a length scale of curvature we may consider 
P = K~ l l A , where K — R^ va£ R^ivae represents the Kretschmann scalar evaluated by contracting Riemann curvature 
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tensor. We must numerically check whether the assumption holds along a given trajectory. 

Binary valued matrix R (/ represents the RP which we get by assigning a black dot where R, 7 = 1 and leaving a white 
dot where Ry = 0. Both axis represent the time period over which the data set (phase space vector) is being examined. 
RP is thus symmetric and the main diagonal is always occupied by the line of identity (LOI). 

Structures in the RP encode surprisingly large amount of information about the dynamics of the system [11]. For 
the purpose of deciding whether the particular trajectory is regular or chaotic the determining factor is a presence 
and amount of diagonal structures in the RP. Diagonal lines in the RP reflect the time segments of the phase space 
trajectory where the system evolves similarly. It captures the epoch when the trajectory runs almost parallel to its 
previous segment i.e. runs inside the e-tube around this segment [5]. Hence integrable systems with regular dynamics 



result in strongly diagonally oriented RPs. On the other hand when the motion is chaotic the diagonal lines are rather 
short (because the trajectory tends to diverge quickly) and more complicated structures are found in the RP. 

The recurrence quantification analysis (RQA) [5] takes number of statistic measures of the recurrence matrix R (; -. 
First of all we define the recurrence rate RR as a density of points in the RP: 

ffl ( e )4l R y(4 (8) 

Then we turn our attention to the diagonal segments in the RP whose length basically draws distinction between 
regularity and chaos. Histogram P(s,l) recording the number of diagonal lines of the length I is formally given as 
follows: 

P(e,l) = £ (l-RM, H (£))(l-R i+1 , j+ ,(£))nR i+ w(4 (9) 

i,j= 1 k=0 

Histogram P(s,l) is used to define the determinism DET as a fraction of the recurrence points which form diagonal 
lines of length at least Z m j n to all recurrence points: 



DET= ^ (10) 



average length of the diagonal line L (where only lines of length at least Z min count): 



and divergence DIV as an inverse of the length of the longest diagonal line L max : 

DIV = — . (12) 

^max 

Since DIV is in its very nature closely related to the divergent features of the phase space trajectory it was originally [8] 
claimed to be directly related to the largest positive Lyapunov characteristic exponent X max - Nevertheless theoretical 
considerations justify the use of DIV as an estimator only for the lower limit of the sum of the positive Lyapunov 
exponents [5]. On the other hand strong correlation between DIV and Amax still arises in numerical experiments [12]. 

It appears useful to perform the analogous statistics also for the vertical (horizontal respectively, since the RP is 
symmetric with respect to LOI) segments in RPs which are generally connected with periods in which the system is 
slowly evolving (laminar states). To this end histogram P(e,v) recording the number of vertical lines of length v is 
constructed as follows: 

P(e,v) = £ (1 -Ry(e))(l -Ry +v (e)) V nR ( - J - + ,(e). (13) 

i,j=l k=0 

In analogy with the diagonal statistics histogram P(e, v) is used to define vertical RQA measures. Laminarity LAM 
is defined as a fraction of recurrence points which form vertical lines of length at least v m ; n to all recurrence points: 

£*_ v . vP(e,v) 
Iv=ivP(e,v) 

trapping time TT is an average length of the vertical line (where only lines of the length at least v m ; n count): 



and finally also the length of the longest vertical line V max might be considered as an analogue of the diagonal L, 



TEST PARTICLE MOTION 



Equations of motion 

Generalized Hamiltonian ("super Hamiltonian") which characterizes the dynamics of the test particle of charge q 
and mass m is given as follows [6]: 

J^=^ v (^-qA tl )(7ty-qA v ), (16) 



where 7fy is the generalized (canonical) momentum and denotes the vector potential related to the electromagnetic 
field tensor as F^ v =A V ^ — A^ jV . 
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FIGURE 1. Classification of possible topologies of off-equatorial potential lobes appearing above the event horizon (bold line) 
of the Kerr black hole immersed into the Wald's test field. See [2] for details. 

The Hamiltonian equations of motion are given in a standard way: 

Ax^/dX = dJtr/dfy, d%/dA = -dJtf/dx 1 *, (17) 

where A = z/m denotes the affine parameter and % the proper time. 

In the case of stationary and axially symmetric systems we immediately identify two constants of motion from the 
second Hamiltonian equation, namely energy E and angular momentum L: 

Jtt=Pt+qA t = —E, n,j ) =p < j ) +qA^=L. (18) 



Effective potential 

We depart from the normalization of the four-momentum 

?"Vv=rW-^)(^-^v) = -m\ (19) 

which is conserved along the Hamiltonian flow in autonomous systems. In stationary and axially symmetric situation 
it allows to express the condition for simultaneous turning point in r, 9 coordinates, i.e. simultaneous zero points of 
u r = p'/m = g rr n r /m and u e = p e jm = g ee %e/m. These points connect in curves which represent natural boundary 
for the test particle motion at given energetic level. This two-dimensional (i.e. related to the motion in two coordinates) 
effective potential describing the motion of the charged test particle takes the following form: 



V eff (r, 0; a, qQ, qB Q , L) = £ + ^ ^ (20) 

denoting 

a = -g", (21) 



p=2[^{L-qAt)-g ,t qA t }, (22) 
7 = -g^ (I - qA^ f - g"q 2 Aj + 2g'*qA, {L — qA§) — 1, (23) 

where we introduce "specific" quantities L = L/m and q = q/m. Setting the specific energy E = E/m of the given 
test particle we obtain isocontours forming the boundary of allowed region which is accessible to its trajectory. Since 
the quantities q, Q and Bo appear only in terms qQ and qBo we only need to specify values of these two products to 
uniquely determine the test particle (provided that remaining parameters are already set). 
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FIGURE 2. The test particle (L = 6M, qBo = M _1 and qQ = 1) is launched from the locus of the off-equatorial potential minima 
r(0) = 3.68 M, 6(0) = 1.18 with u r (0) = and various values of the energy E. In the left panel we set E = 1.58 and we observe 
ordered off-equatorial motion. For the energy of E = 1 .65 cross-equatorial regular motion is observed (middle panel). The trajectory 
occupies only a part of allowed potential lobe, regardless the length of the integration period. Finally in the right panel with E = 1 .75 
we observe irregular motion whose trajectory would ergodically fill whole allowed region after the sufficiently long integration time. 
We show that the motion is chaotic in this case. Spin of the black hole is a = 0.9 M and its event horizon is depicted by the bold 
line. 

In order to locate the halo orbit and related off-equatorial potential lobe we search for the local minima of the 
potential V e ff. Direct approach based on the investigation of the conditions = 0, -4# = and the determinant of 
the Hessian matrix being positive, leads to the set of very intricate algebraic equations. However, when we reformulate 
these conditions using the force formalism we obtain cubic equation which we are able to solve in general. See [3] 
for details. We found [2] that potential minima may appear above the event horizon of the black hole, i.e. that off- 
equatorial stable orbits are allowed in this setup - see Fig. 1. Besides the off-equatorial lobes the potential may form 
another remarkable structure - endless potential valley of almost constant depth which runs parallelly to the symmetry 
axis. Thus the particle may escape to infinity from the equatorial plane if its parameters are suitably chosen. Such orbits 
may form highly collimated jet-like structure and may also provide a charge separation mechanism as was recently 
reported in [13]. Motion in the potential valleys we further discuss in [1]. Here we concern ourselves with the motion 
in the closed potential lobes only. 

REGULAR AND CHAOTIC MOTION IN OFF-EQUATORIAL LOBES 

In this section we shall apply above described recurrence analysis to the chosen set of orbits. From our classification 
of the off-equatorial potential lobes [2] which may appear above the horizon of the Kerr black hole immersed in the 
asymptotically uniform magnetic field (setup described by equations (1), (3) and (4)) we choose the type Id from 
Fig. 1 for our analysis. In this case the off-equatorial local minima of the effective potential (20) are surrounded by the 
potential lobes which merge through the equatorial plane once the energy level of the equatorial saddle point is reached 
and form closed lobe which extends symmetrically above and below the equatorial plane. Increasing the energy even 
more we would eventually reach the level when this lobe breaks symmetrically in its bottom and upper parts allowing, 
in principle, the particle to escape to the infinity in the axial direction. Further increase in the energy leads to the 
opening of the lobe towards the event horizon through the saddle point in the equatorial plane. Described situation is 
captured in Fig. 2. 
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FIGURE 3. For the sake of comparison we present the surface of section and recurrence plot for a regular trajectory in the fully 
integrable system of charged test particle (E = 0.9965, L = 6M, q = 10 17 , r(0) = 70 A/, 9(0) = e sec tion = w/2 and u r (0) = 0) in 
the pure Kerr-Newman spacetime (Q = 5 x 10~ 18 , a = 0.9 M) endowed with the fourth Carter's constant of motion «Sf [6]. Long 
diagonals parallel to the LOI are a general hallmark of regularity in the RPs. 
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FIGURE 4. Regular off-equatorial motion of the charged test particle (E = 1.578, L = 6M, qQ = 1, qB Q = AT 1 , r(0) = 3.68M, 
m''(0) = and 0(0) = ^section = 1.18) on the Kerr background (a = 0.9 M) with the Wald's test field. In the Poincare surface of 
section we distinguish u e > (black point) from u e < (grey point). 



Numerical integration of the Hamilton's equations (17) is carried out using the multistep Adams-Bashforth-Moulton 
solver of variable order. In several cases when higher accuracy is demanded we employ 7-8th order Dormand-Prince 
method. Initial values of non-constant components of the canonical momentum 7t r (0) and 7tg(0) are obtained from 
u r (0) (which we set) and u 6 (0) which is calculated from the normalization condition g^u^Uy = - 1 where we always 
choose the non-negative root as a value of u e (0). CRP ToolBox 1 is used for the construction of RPs and evaluation of 
RQA measures. 

In Fig. 3 we show how does the regular trajectory in the fully integrable system of the charged test particle in the 
Kerr-Newman spacetime appear in the Poincare surface of section and in the recurrence plot. We observe that in this 
case the RP exhibits purely diagonal features as expected. 



1 CRP ToolBox for Matlab: http://www.agnld.uni-potsdam.der marwan/toolbox/ 
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FIGURE 5. Regular cross-equatorial motion of the charged test particle exposed to the Wald's field which only differs from the 
previous case of Fig. 4 by increasing energy to E = 1.65 for which both off-equatorial potential lobes are already merged via the 
equatorial plane. 




FIGURE 6. Chaotic motion of the test particle which appears when the energy increases to E = 1 .75 (while all other parameters 
are kept at values of Figs. 4 and 5). 



Then in Fig. 4 we turn attention to our non-integrable system of the charged test particle being exposed to the Wald's 
field on the Kerr background. At given energy level of E = 1 .578 the particle remains trapped in the off-equatorial lobe. 
Its motion is regular and diagonal structures in the RP typical for the trajectories in integrable systems are preserved, 
though the pattern is slightly different. 

Increasing the energy to E = 1 .65 we obtain the lobe already merged via the equatorial plane. To some surprise 
the merging itself which occurs at E w 1 .59 is not reflected by the change of the dynamic regime - although the test 
particle can newly cross the equatorial plane, its motion remains regular (see Fig. 5). Nevertheless the RP proves its 
dynamics to be more complex than that of the integrable system of Fig. 3. Diagonal structures are preserved, though 
the pattern is more complicated. 

However, at E = 1 .75 we obtain typical chaotic motion Fig. 6. In the surface of section we observe that the trajectory 
fills all allowed region (hypersurface given by the values of integrals of motion). Diagonal lines in the RP are seriously 
disrupted and complex large scale structures appear which are characteristic indications of deterministic chaos. 

We conclude that the motion in purely off-equatorial lobes (i.e. before merging of the lobes via equatorial plane) 
is regular. Surprisingly the regularity of motion is not violated when the lobes merge and volume of the accessible 




portion of phase space suddenly doubles. However, increasing further the energy eventually brings the system to the 
chaotic regime of motion. This route to chaos was documented by both the standard Poincare surfaces of section and 
the recurrence plots. RPs prove to act as an alternative tool to surfaces of section as they allow to distinguish between 
the regular and the chaotic regime of motion by apparent changes in their patterns. 

So far we know that in our set of orbits the transition between the regular and the chaotic regime of the motion 
occurs somewhere in the range E 6 (1.65,1.75). Qualitative methods based on the visual survey of the Poincare 
surfaces of section and the recurrence plots do not allow to localise this transition. To this end we employ recurrence 
quantification analysis (RQA) and observe the behaviour of RQA measures (definitions (8)-(15)). We calculate 200 
trajectories with energies equidistantly spread over the range (1.65, 1.75) while other parameters of the system are 
fixed at values used in Figs. 4-6. In Figs. 7-10 we observe that all queried RQA measures exhibit sudden change in its 
behaviour at E « 1 .685. This is where the dynamic transition between the regimes occurs. Transition from the regular 
motion to the chaos is thus detected not only by diagonal RQA measures but also by the vertical ones. 
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FIGURE 9. RQA measures derived from vertical lines LAM and TT also clearly signalize the onset of chaos at E ~ 1 .685. 
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FIGURE 10. Length of the longest vertical line V ma x generally rises and begins to fluctuate strongly as the chaos sets on at 
1.685. 

CONCLUSIONS 

In this contribution we demonstrated that the recurrence plots and the recurrence quantification analysis are simple, yet 
powerful tools which allow one not only to decide whether the dynamic regime of motion is regular or chaotic but also 
to locate (in terms of some control parameter - energy E in our case) the transition between these regimes with a good 
precision. General relativistic context highlights the advantages of the recurrence analysis over standard methods as 
it operates "locally" (on the small length scale of e) which allows for various profound computational simplifications 
when evaluating distances in the phase space. Complex treatment which involves evaluation of the geodesic line and 
measuring its length may thus be avoided. Major drawback (cost we pay for its simplicity) of RPs and RQA is the lack 
of invariance (dependence on the threshold e, / m i n , v m j n , choice of the norm etc.). We conclude that RPs themselves 
may act as an alternative tool to Poincare surfaces of section and RQA measures are able to detect the transitions 
between the dynamic regimes. 
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